Optical excitations of Si by time-dependent density-functional theory based on the 

exact-exchange Kohn-Sham band structure 
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We calculate the imaginary part of the frequency-dependent dielectric function of bulk silicon by 
applying time-dependent density-functional theory based on the exact-exchange (EXX) Kohn-Sham 
(KS) band structure and the adiabatic local-density approximation (ALDA) kernel. The position 
of the Ei absorption peak calculated with the EXX band structure at the independent-particle 
level is in excellent agreement with experiments, which demonstrates the good quality of EXX 'KS 
quasiparticles'. The excitonic E\ peak that is missing at the independent-particle level remains 
absent if two-particle interaction effects are taken into account within the time-dependent LDA, 
demonstrating the incapability of the ALDA kernel to describe excitonic effects. 
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I. INTRODUCTION 

While the Kohn-Sham (KS) density-functional theory 
(DFT) 0, ||, ||] has become the dominant method of first- 
principles investigations of energetic properties of solids, 
its scope in the treatment of electronic excitations has 
been rather limited due to deficiencies of the conventional 
local density approximation (LDA) and its semilocal ex- 
tensions. One well-known case is the band gaps calcu- 
lated within the LDA which are substantially smaller 
than experimental values 0). Significant advancement 
has been made during the past decades in devising im- 
proved exchange-correlation functionals starting from the 
LDA J3j, but these semilocal approximations still fail 
to correct the so-called 'band-gap problem'. Thus for 
the accurate treatment of band structures one at present 
usually leaves the domain of DFT and employs a many- 
body method based on the GW approximation (GWA) 
to quasiparticle energies ||, |j| . 

For the case of optical excitations in which bound ex- 
citon states can play a significant role, it is further re- 
quired to go beyond effective one-particle approaches. In 
principle, this can be done within the time-dependent 
(TD) DFT §, § or by solving the Bethe-Salpeter equa- 
tion (BSE) for two-body Green's functions M. How- 
ever, the lack of a proper exchange-correlation functional 
again remains as the hindrance to the TDDFT route, and 
the GWA-BSE currently represents the state-of-the-art 
method for such purpose [H], [ll[] . 

In this regard recent numerical realization of the exact- 
exchange (EXX) KS method for solids Jli| and 
molecules fi| g H || |§ provides an interesting 
opportunity: According to theoretical and numerical ev- 
idence concerning the physical meaningfulness of the KS 
potential and KS eigenvalues [^l], |2^, p3| , realistic band 
structures for semiconductors and eigenvalue spectra of 
molecules obtained within the EXX suggests a promis- 
ing basis for the study of electronic excited states within 
DFT. In this work, we consider this point further by cal- 
culating the absorption spectrum of bulk silicon within 



TDDFT based on the EXX band structure and the adi- 
abatic LDA (ALDA) exchange-correlation kernel. Sili- 
con here serves as a representative semiconductor system, 
whose first absorption peak, the E\ peak, is known to be 
strongly excitonic in its character. The large excitonic 
contributions to the E\ peak has been proposed theoreti- 
cally and experimentally, and a semiempirical calculation 
confirming this point has been performed two decades 
ago ||. However, truly first-principles calculations ap- 
peared only recently based on the GWA-BSE || [n). 
A different approach to investigate optical spectra within 
a DFT framework was recently explored in Ref. |24|, where 
the domain of conventional DFT was left and time- 
dependent current DFT was invoked. As a first step to- 
wards an accurate computation of optical spectra within 
TDDFT, we will examine the performance of the EXX- 
TDLDA approach in this paper. 



II. THEORY 

To obtain the absorption spectrum of solids, we com- 
pute the imaginary part of the frequency-dependent 
macroscopic dielectric function t(uj) = t\{uj) + ie2(w) de- 
fined as l2l |2t| 



e(uj) = lime 1 (q + G, q + G'; uj)\ 

q — »q 1 — u 



where 

e-^q + G, q + G'; u) = 5 G ,G> + F H (q + G) 



x X (q + G,q+G';^), 



(1) 



(2) 



with Fff(q+G) = 47r/|q-|-G| 2 (Hartree atomic units are 
used throughout the paper. The symbols q + G refer to 
plane waves with wave vectors q + G, e l ( q+G )' r , i.e., all 
formulas are given in the plane wave representation) . The 
linear response function x that describes the full response 
of the first-order density change Sn for the given external 
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perturbation 5v ex t, 

<5n(q+G;w) = x(q + G, q + G'; (5u e xt(q+ G';w), 

G' 

(3) 

can be computed either via the BSE or TDDFT. Within 
TDDFT the linear response function x is given in matrix 
representation as |], [j 

r i _1 

x(q; w) = [l - xo(q; ui){F H (q) + F sc (q; w)}J xo(q; w). 

(4) 

where xo is the noninteracting KS linear response matrix. 
It is given in terms of the occupied and unoccupied KS 
orbitals <^>„k and </> c k and their eigenvalues e„k and e c k 
through 



Xo(q + G,q + G';^) = — 



tensor 



E 



(v, k\e-^ +G ^ r \c, k + q) (c, k + q|e*( q+G ')- r ' \v, k) 
e-uk - e c (k+ q ) + iS 



(c, k|e- l ( q+G )- r |w, k + q)(v, k + q\e l ^ +G '> r ' |c, k) 



e uk e c(k+q) 



id 



(5) 

where Nil is the crystal volume. In Eq. (|J) F^c denotes 
the exchange-correlation kernel, the frequency-dependent 
functional derivative of the exchange-correlation poten- 
tial. 

Thus, the macroscopic dielectric function e can be com- 
puted by first calculating e _1 in the reciprocal-space us- 
ing Eqs. (||) and (||), and next taking the inverse of the 
q -> limit of its G = G' = element [Eq. (Q)]. How- 
ever, the ALDA kernel F A C LDA = S 2 E^ A /(6n Sn') is 
local in real space and frequency- independent, which re- 
sults in a reciprocal-space representation for F xc that is 
independent of q and to, F A c LDA (q + G,q + G';w) = 
F xc LDA ( G ~ G )■ This incorrect q -> behavior of 
the ALDA kernel then simplifies the calculation of e(w): 
Writing x, Fjj, and F xc (from now on the superscript 
"ALDA" on the kernel is suppressed for notational sim- 
plicity) for q — > in block form in terms of the "head" 
[G = G' = (00)], "wings" [G = 0; G' ^ (01) and 
G ^ 0;G' = (10)], and the "body" [G ^ 0; G' ^ 
(11)] as 



xo(g;w) = Q 



F H (q) = Q- 



x°°( w ) xo» 
xJV) xo»' 1 



4tt 








Q~ 



(0) 



where Q 00 = g; Q 01 = g 10 = ; g n = 6 G , G ', and 
performing matrix multiplications explicitly (with q- 
singularities treated properly at each stage) yields a 3 x 3 



e( W ) = 1 - 47r X ^(c) - 47r X g 1 M(^ 1 + 



(7) 



with which matrix multiplications including head and 
wings with explicit q-dependence are avoided. Eq. (j?]) 
simply represents the independent-particle-level absorp- 
tion spectrum if Fj 1 and F^ are ignored altogether, 
while it is reverted to the so-called RPA dielectric func- 
tion formula if only F^ is neglected. For cubic systems, 
such as Si considered here, the macroscopic dielectric ten- 
sor is diagonal. 



III. COMPUTATIONAL METHODS 

For the calculation of the dielectric function, one needs 
wavefunctions (j)^ and eigenvalues e^k at the chosen set 
of k-points [See Eq. (||)]. To generate these, we first 
performed the LDA and EXX self-consistent ground- 
state calculations at the experimental lattice constant, 
a = 5.43A for Si, and obtained the KS potential Vks 
which serves as the only input to the TDDFT calcula- 
tion of optical spectra. The EXX calculation was car- 
ried out using ten special-k points, orbital kinetic energy 
cutoffs of 12.5 Ha, and nonlocal pseudopotentials in the 
separable form of Kleinman-Bylander j2?j] . We refer the 
reader to Ref. [l4] for further details of EXX calculations. 
In addition, because the ALDA kernel F A C LDA depends 
only on the density, we constructed it at the end of the 
ground-state calculations when fully converged densities 
were available. 

With the converged wff x , v^g A and F A C LDA in hand, 
we moved on to the computation of the dielectric function 
and first solved the KS equations once again at a much 
larger number of k-points but with the same energy cutoff 
12.5 Ha to generate xo- The calculation of the head and 
wings of xo needs special care since the q — > limit of 
(i, k|e _lq ' r |j, k+q) should be handled properly. Applying 
per turbation theory one obtains the following equation 



lim(i, kle 



|i,k + q) =q 



,k|p + [V NL ,ir]\j,k) 



(8) 



which implies that, due to the extra powers of the energy 
denominator that appear in the head Xq (q;w)> compu- 
tation parameters for the convergence are more or less 
determined by those necessary for the head element. We 
use the kinetic energy cutoff 7 Ha and 10 conduction 
bands for xo which was sufficient to reach the conver- 
gence of £2 (w). 

For the k-point sampling, we chose to adopt a 
regular uniform grid, with which one could evaluate 
limg^o^; k|e _lq r |j, k + q) by a finite difference approach 
without resorting to the analytic replacement of Eq. (0). 
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In this case one had to pay the price of performing three 
additional diagonalizations at k + qe a (a = x,y,z) for 
each k-point. Therefore, in this paper we report results 
obtained by treating the q — » limit according to Eq. (||) 
with ignoring the nonlocal pseudopotential contributions 
due to the commutator [VivL,*r] term. In test calcu- 
lations with a small number of k-points we found that 
the neglect of the nonlocal pseudopotential contributions 
has little effect on the positions of absorption peaks while 
the peak heights are somewhat reduced. We slightly shift 
the uniform k-mesh in order to break any present sym- 
metry thus achieve a more effect k-point sampling [p_ 1| . 
After checking the convergence behavior by increasing 
the number of k-points and changing the shifting vector, 
we used a shifted 17x17x17 regular grid which gives 
the convergence of £2 up to the frequency range of at 
least 6 eV. The use of a shifted uniform k-mesh results 
in slightly non-identical diagonal elements of the dielec- 
tric tensor whose differences have been also used to check 
the convergence with respect to the number of k-points. 



IV. RESULTS AND DISCUSSIONS 

We first consider the absorption spectrum calculated 
at the LDA and EXX independent particle levels. In 
Fig. LDA and EXX absorption spectra are shown to- 
gether with experimental jp8| and GWA || [l(| curves. 
Note that the two GWA results are rather different due 
to operational ambiguities resulting from the post-DFT 
nature of the GWA g, while the LDA and EXX which 
strictly reside in the KS DFT scheme do not suffer from 
such problems. Within the LDA, the absorption spectra 
is shifted to the lower frequency region by about 1 eV 
due to its well-known underestimation of bands gaps. In 
addition, the E\ peak height is significantly underesti- 
mated, while that of the E 2 peak is much overestimated. 
In contrast, the EXX absorption spectrum shows a much 
better overall agreement reflecting realistic descriptions 
of the band structure within the EXX |14J: Positions 
of the absorption edge and the E 2 peak are in excellent 
agreement with the experiment, and the E% peak posi- 
tion is even slightly better than the GWA one. Also the 
height of the E2 peak is smaller than the LDA value and 
is much closer to the experimental data. One noticeable 
disagreement is the absence of the E\ peak as in the LDA 
spectrum. This, however, does not reflect the deficiency 
of the EXX approach, but the excitonic nature of the E\ 
peak H |l(| [ll| which cannot be described at the inde- 
pendent particle level as can be also seen in the GWA 
curves. 

We feel that the excellent features of the EXX absorp- 
tion spectrum warrant more discussions. Perturbation 
theory along the adiabatic connection shows that KS 
eigenvalue differences are well-defined approximations of 
excitation energies [^IJ. Indeed Umrigar et al. recently 
found for several atoms that quasi-exact KS eigenvalue 
differences are in surprisingly good agreement with ex- 



perimental excitations J22[. In addition, applying the 
quasi-exact KS potential to the calculation of the quan- 
tum defects for the Ne atom, Al-Sharif et al. argued 
that the KS potential of a neutral atom (with TV elec- 
trons) corresponds to the Green-function "potential" of 
its daughter cation (TV — 1 electrons), or equivalently the 
KS equation for TV electrons corresponds to the Dyson 
equation where the reference ground state is chosen with 
TV — 1 electrons p3| |. The appealing description of exci- 
tations by KS eigenvalue differences, in contrast to the 
treatment of ionizations and electron affinities in the 
many-body (GWA) approach, has not drawn much at- 
tention due to the lack of realistic exchange-correlation 
approximations which are suitable for practical calcu- 
lations. Exchange-correlation potentials derived from 
the LDA and other semilocal exchange-correlation en- 
ergy functionals do not cancel Coulomb self-interactions, 
decay exponentially rather than as —1/r for localized sys- 
tems, and inherently lack the term that results from the 
functional derivative of the pair-correlation function p9[ . 
The self-interaction-free EXX method corrects such prob- 
lems and provides accurate KS potentials and eigenval- 
ues. Then, we propose that the good absorption spec- 
trum from such EXX method represents further evi- 
dence for the TV-electron excitations picture within the 
KS scheme. 

Now we consider the effects of the inclusion 
of two-particle interactions via TDLDA on top 
of the independent-particle EXX result. The 
Hartree-only or RPA (EXX+TDH) and full TDLDA 
(EXX+TDH+ALDA) spectra are shown together with 
the EXX independent-particle spectrum in Fig. |^. 
Granted the overall good feature of the EXX spectrum, 
the most interesting question is that whether the ALDA 
can generate the excitonic E\ peak. Fig. || shows 
that this is not the case. TDH alone suppresses the 
EXX spectrum and shifts the E2 peak to a slightly 
higher frequency region, while adding the ALDA coun- 
teracts such effects and more or less restores the EXX 
independent-particle spectrum. In Fig. || we also display 
the TDLDA absorption spectrum based on the LDA band 
structure (LDA+TDH+ALDA) to show that the LDA 
independent-particle result (see Fig. Q) is again insignif- 
icantly changed by TDLDA effects. 



V. CONCLUSIONS 

In this work, by calculating the frequency-dependent 
dielectric function of bulk Si we studied the quality of 
EXX one-particle band structure and the description of 
excitonic effects within the TDLDA. We showed that, at 
the EXX independent-particle level, the position of the 
Ei peak is in good agreement with experimental data 
and is even slightly better than that obtained in the 
GWA approach. We argued that this is because the EXX 
approach, unlike procedures employing the LDA and 
semilocal density-functionals, yields realistic KS orbitals 
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and eigenvalues which accurately model electronic exci- 
tations at a fixed particle number. Employing TDDFT 
within the ALDA on top of the EXX band structure how- 
ever results in only a minor modification of the EXX 
spectrum and fails to produce the excitonic E\ peak. 
Given the fact that proper inclusion of excitonic effects 
is mandatory for the complete description of optical ex- 
citations, devising a TDDFT scheme that overcomes the 
deficiency of the ALDA kernel will be an interesting task. 
Work along this line using the nonadiabatic and nonlocal 
EXX kernel is in progress (pOf. 
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FIG. 1: Absorption spectrum of bulk silicon at the one-particle LDA and EXX levels. Experimental curve is from Ref. 
Two GWA results (GWA1 from Ref. [To] and GWA2 from Ref. |^) are also shown for comparison. 
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FIG. 2: Absorption spectrum of bulk silicon at the EXX, EXX+TDH, EXX+TDH+ALDA, and LDA+TDH+ALDA levels. 
Experimental data are from Ref. EM 



